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Abstract 

In  this  paper  we  develop  a  mathematical  model  for  the  rapid  production  of  large 
quantities  of  therapeutic  and/or  preventative  countermeasures.  We  couple  equations 
for  biomass  production  with  those  for  vaccine  production  in  shrimp  that  have  been 
infected  with  a  recombinant  viral  vector  expressing  a  foreign  antigen.  The  model 
system  entails  both  size  and  class-age  structure. 
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1  Introduction 


We  consider  a  novel  approach  for  developing  a  stable  operational  platform  for  the  rapid 
production  of  large  quantities  of  therapeutic  and/or  preventative  countermeasures. 
The  ideas  developed  here  can  also  serve  as  the  foundations  in  designing  an  economical 
platform  for  the  production  of  complex  protein  therapeutics  to  replace  mammalian 
cell  culture  production  methods  used  in  the  pharmaceutical  industry.  This  approach 
involves  recruiting  the  biochemical  machinery  in  an  existing  biomass  for  the  production 
of  a  vaccine  or  antibody  by  infection  using  a  virus  carrying  a  passenger  gene  for  the 
desired  countermeasure.  While  our  motivation  derives  from  efforts  related  to  first 
response  to  deliberate  bio  toxic  attacks  on  populations,  the  models  we  develop  may  also 
have  use  in  designing  prophylactic  production  systems  against  epidemics  originating 
naturally  in  populations  which,  without  intervention,  might  result  in  pandemics.  While 
our  model  is  specific  to  virus  growth  and  vaccine  production  in  shrimp,  the  implications 
for  other  crustaceans  are  obvious.  And  of  course  the  shrimp  models  we  investigate  can 
serve  as  a  foundation  for  understanding  viral  progression  in  other  species  important  to 
marine  agriculture. 

Our  goal  is  to  model  a  system  wherein  one  uses  shrimp  as  a  scaffold  organism  to 
produce  biological  countermeasures.  In  such  a  system  one  might  first  stock  shrimp 
postlarvae  and  allow  them  to  grow  normally  in  a  controlled  environment.  Then  one 
infects  them  with  a  recombinant  viral  vector  (e.g.,  recombinant  Taura  Syndrome  virus 
or  rTSV)  expressing  a  foreign  antigen,  resulting  in  vaccine  production  in  live  infected 
shrimp. 

To  mathematically  demonstrate  the  feasibility  of  this  approach  we  consider  a  hybrid 
model  of  the  shrimp  biomass/countermeasure  production  system  which  has  two  com¬ 
ponents:  biomass  production,  and  production  of  countermeasure  (antibody/vaccine). 
We  feed  the  output  of  the  biomass  production  model  as  input  to  the  vaccine  production 
model.  For  initial  investigations  the  amount  of  vaccine  produced  is  assumed  equal  to 
the  total  infected  biomass.  Thus,  the  vaccine  production  model  will  essentially  follow 
the  course  of  the  viral  dynamics  in  shrimp. 

There  is  considerable  literature  on  shrimp  growth  dynamics  (e.g.,  [10,  11,  17,  24]) 
where  mathematical  descriptions  of  growth  curves  are  fitted  to  data  obtained  for  a 
specific  subspecies  of  shrimp  and  no  structural  information  is  utilized.  However,  there 
is  little  information  in  the  literature  on  modeling  the  dynamics  of  shrimp  at  the  pop¬ 
ulation  level.  In  [9],  the  authors  incorporate  abiotic  factors  such  as  temperature  into 
a  logistic  type  growth  equation  to  model  biomass  dynamics  in  shrimp;  in  such  models 
it  is  assumed  that  all  individuals  are  identical  in  characteristics  and  behavior.  How¬ 
ever,  disregarding  structure  in  constructing  mathematical  models  for  the  dynamics  of 
shrimp  is  unrealistic,  since  shrimp  have  size  dependent  characteristics  and  responses  to 
external  environment.  The  model  in  this  paper  is  based  on  the  classical  McKendrick- 
von-Foerster/Sinko-Streifer  size-structured  population  equations  [18,  21,  23]  with  mass 
as  the  structure  variable,  i.e.,  we  equate  the  size  variable  with  the  mass  in  our  model. 

There  also  appears  to  be  a  dearth  of  literature  on  modeling  epidemics  in  shrimp 
populations.  In  [20]  the  authors  develop  a  non-structured  five  compartment  epidemic 
model  of  TSV  that  includes  a  Reed-Frost  transmission  process  in  closed  populations 
of  shrimp  ( Litopenaeus  vannamei).  However,  as  in  the  case  of  the  biomass  model, 
structure  can  play  an  important  role  in  the  study  of  viral  epidemiology  in  shrimp. 
Moreover,  experimental  results  [14]  suggest  that  the  mortality  rate  in  acutely  infected 
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shrimp  depends  on  the  length  of  time  that  the  shrimp  remain  acute.  Also,  individuals 
in  the  latent  phase  have  varying  residency  times  before  they  progress  into  the  acute 
phase.  To  incorporate  all  of  these  features,  we  attempted  to  model  the  progression  of 
TSV  in  shrimp  in  a  system  of  delay  PDEs.  However,  it  is  difficult  to  correctly  account 
for  the  different  residency  periods  of  individual  shrimp  in  this  fashion  as  the  size  of  the 
shrimp  is  a  function  of  time.  Instead  of  tracing  back  in  time  to  incorporate  delays,  a 
different  approach  involves  recording  the  variable  residency  times  in  the  different  stages 
by  introducing  a  new  variable  which  we  call  the  class  age  of  an  individual.  The  class  age 
of  an  individual  in  a  given  stage  represents  the  length  of  time  that  the  individual  spends 
in  that  stage.  Ours  is  the  first  attempt  at  a  mathematical  model  that  incorporates  size 
and  class  age  to  study  the  progression  of  TSV  infection  in  shrimp.  A  similar  approach 
has  been  used  previously  to  investigate  a  linear  cell  population  model.  In  such  models, 
cells  are  assumed  capable  of  simultaneous  proliferation  and  maturation  where  in  the 
proliferating  phase,  cells  are  committed  to  undergo  cell  division  some  time  units  after 
entering  this  phase  [1] .  There  are  several  papers  that  investigate  structured  population 
models  with  multiple  internal  variables,  for  example,  see  [15,  22]. 

The  outline  of  the  paper  is  as  follows:  In  Section  2  we  plot  the  growth  data  observed 
in  several  different  shrimp  populations  cultured  in  super-intensively  stocked  raceways 
production  systems  at  the  Waddell  Mariculture  Center  in  Bluffton,  SC.  We  use  regres¬ 
sion  analysis  to  fit  linear,  cubic  and  exponential  growth  rate  models  to  this  data.  In 
Section  3  we  formulate  a  biomass  production  model  for  the  normal  growth  of  shrimp 
and  present  results  of  numerical  simulations  of  this  model,  while  in  Section  4  we  con¬ 
struct  a  three  compartment  vaccine  production  model  and  analyze  this  model  using 
the  method  of  characteristics.  We  also  construct  a  discrete  scheme  for  the  vaccine 
production  model.  Numerical  simulations  for  the  coupled  biomass-vaccine  production 
model  are  discussed  in  Section  5. 

Successful  implementation  of  this  shrimp-based  expression  system  could  potentially 
provide  an  inexpensive  alternative  method  to  mammalian  cell  culture  for  making  com¬ 
plex  recombinant  proteins.  In  this  context,  complex  recombinant  proteins  are  ones 
that  undergo  post-translational  modifications  such  as  glycosylation  that  cannot  easily 
be  duplicated  in  bacterial  and  yeast  protein  expression  systems.  The  shrimp  expression 
system  could  dramatically  lower  the  costs  of  production  of  these  complex  proteins. 


2  Data  Fitting 

We  have  fitted  linear,  cubic  and  exponential  solutions  x  for  individual  growth  rate 
dx 

models  —  =  g(x,  t )  (corresponding  to  g  constant,  quadratic  and  exponential)  to  the 
data  for  the  mean  weight  of  the  shrimp  measured  during  a  series  of  grow  out  runs  in  a 
super-intensively  stocked  demonstration  commercial  scale  raceway  production  system 
[7]  at  the  Waddell  Mariculture  Center.  Growth  data  is  based  on  population  sampling 
from  a  series  of  demonstration  runs  in  a  new  advanced  greenhouse  enclosed  biosecure 
system  developed  for  the  culture  of  marine  shrimp  Litopenaeus  vannamei.  The  Waddell 
system  provides  for  the  culture  of  shrimp  at  very  high  densities  without  water  exchange. 
A  dense  microbial  floe  develops  in  the  water  column  providing  for  the  recycling  of  waste 
material  within  the  system  to  maintain  water  quality  and  enhance  shrimp  growth. 
By  applying  filtration,  oxygen  injection  and  well  designed  feeds  and  feeding  regimes, 
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genetically  selected  specific  pathogen  free  strains  of  shrimp  can  be  cultured  in  this 
system  at  high  densities  with  good  survival  and  excellent  growth. 

The  various  raceway  runs  initially  had  different  stocking  densities.  The  mean  weight 
x  (in  gm)  of  the  shrimp  were  recorded  on  different  days  over  several  weeks.  The  results 
of  the  regression  analysis  are  presented  below. 

In  raceway  runs  6  and  9,  the  mean  stocking  weight  is  1  gm,  and  linear  fits  give  fair 
approximations  to  the  data  obtained  from  these  runs  as  observed  in  Figures  1  and  3, 
respectively.  We  can  also  see  in  these  two  figures  that  a  cubic  fit  gives  a  much  better 
approximation  to  the  data  than  the  corresponding  linear  fits.  Run  10  was  stocked  with 
shrimp  that  had  a  mean  weight  of  4  gm,  and  a  linear  fit  seems  to  be  appropriate  in 
this  case  as  seen  in  Figure  4  (left).  Run  8  was  stocked  with  shrimp  that  had  a  mean 
weight  of  0.01  gm.  From  Figure  2  we  can  see  that  a  linear  fit  does  not  do  justice  to  the 
data.  However,  a  cubic  fit  works  well  in  this  case.  In  all  these  figures,  the  horizontal 
axis  represents  time  t  in  days  and  the  vertical  axis  represents  the  mean  weight  x  of 
the  shrimp.  Each  figure  also  records  the  equation  for  the  corresponding  linear,  cubic 
or  exponential  fit  as  the  case  may  be. 

We  can  conclude  that  a  linear  fit  is  a  good  approximation  to  the  data  in  the  range 
1-20  gm,  but  a  linear  fit  does  not  approximate  the  data  in  the  interval  0-20  gm.  In  [12], 
the  authors  have  observed  that  growth  of  early  postlarvae  was  exponential  rather  than 
linear.  However,  we  were  unable  to  use  an  exponential  function  to  fit  the  entire  data, 
i.e.,  in  the  range  0-20  gm.  For  example,  Figure  4  (right)  depicts  the  best  exponential 
fit  to  the  data  obtained  from  run  8.  We  can  see  that  this  fit  does  not  give  a  good 
approximation  to  the  data.  Instead  a  cubic  fit  seems  to  yield  a  good  approximation  as 
seen  in  Figure  2  (right). 

Using  a  linear  fit  for  the  mean  weight  x  of  the  shrimp  as  a  function  of  time  t 
implies  that  the  growth  function  g(x,  t )  is  a  constant  in  the  interval  1-20  gm.  This 
is  the  approximation  that  we  have  used  to  perform  our  simulations.  Since  a  cubic  fit 
(for  x  as  a  function  of  time  t)  gives  a  far  better  approximation  to  the  data,  as  seen  in 
Figures  1,  2,  and  3,  the  growth  function  g(x,t)  would  then  be  a  quadratic  function  of 
time  t.  We  have  as  yet  not  used  such  an  approximation  in  our  simulations  since  we 
have  assumed  a  mean  stocking  weight  of  1  gm. 


linear  fitting:  raceway  6 


cubic  fitting:  raceway  6 


Figure  1:  Linear  and  cubic  fitting  to  the  data  of  raceway  run  6,  where  the 
mean  stocking  weight  of  the  shrimp  is  1  gm. 
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linear  fitting:  raceway  8  cubic  fitting:  raceway  8 


Figure  2:  Linear  and  cubic  fitting  to  the  data  of  raceway  run  8,  where  the 
mean  stocking  weight  of  the  shrimp  is  O.Of  gm. 


linear  fitting:  raceway  9  cubic  fitting:  raceway  9 


Figure  3:  Linear  and  cubic  fitting  to  the  data  of  raceway  run  9,  where  the 
mean  stocking  weight  of  the  shrimp  is  f  gm. 


linear  fitting:  raceway  10  exponential  fitting:  raceway  8 


Figure  4:  (left)  Linear  fitting  to  the  data  of  raceway  run  fO,  where  the  mean 
stocking  weight  of  the  shrimp  is  4  gm.  (right)  Exponential  fitting  to  the  data 
of  raceway  run  8,  where  the  mean  stocking  weight  of  the  shrimp  is  O.Of  gm. 

3  Biomass  Production  Model 

We  present  a  classical  size-structured  population  model  with  mass  as  the  structure 
variable.  The  factors  considered  here  that  affect  the  total  biomass  are  the  growth  and 
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mortality  rates  of  normal  shrimp.  Since  we  have  a  controlled  environment  for  growth, 
we  may  assume  that  the  growth  and  mortality  rates  of  normal  shrimp  are  not  affected 
by  changes  in  environmental  factors  such  as  temperature,  as  well  as  other  factors  such 
as  population  density. 

Growth  in  shrimp  consists  of  periods  involving  molting  separated  by  internrolt  peri¬ 
ods  where  no  external  growth  occurs  (see  [8]  for  more  details).  This  makes  the  growth 
in  shrimp  a  discontinuous  process.  However,  viewed  over  a  long  period  of  time  we  can 
approximate  the  growth  as  a  continuous  process.  Hence,  we  assume  the  individual 
growth  rate  g(x,  t )  to  be  a  continuous  function  of  mass  x  and  time  t  as  depicted  in 
Section  2.  The  time  interval  for  normal  growth  that  we  consider  here  is  not  sufficiently 
long  for  the  shrimp  postlarve  to  grow  to  adult  shrimp  and  thus  we  do  not  consider 
reproduction  in  the  biomass  and  vaccine  production  models. 

Based  on  the  above  discussion  our  biomass  production  model  is  given  by 

ut  +  (g(x,  t)u)x  +  m(x ,  t)u  =  0,  (x,  t )  £  (0,  xmax]  x  (0,  TB\ , 

n(0,  t)  =  0,  t  £  (0,  TB\,  (3.1) 

u(x,  0)  =  u°(x),  x  £  [0,  xmax]. 

Here  u(x,t )  denotes  density  of  individuals  (number/unit  mass)  having  mass  x  in 
gms  at  time  t.  The  growth  rate  is  denoted  by  g(x,  t )  (mass/unit  time),  and  the  function 
m(x,  t )  denotes  the  mortality  rate  (1/unit  time).  The  initial  population  density  is  given 
by  u°(x)  and  xmax  is  the  maximum  mass  of  shrimp  in  the  time  interval  from  0  to  a 
final  time  TB. 

There  are  numerous  papers  that  consider  wellposedness  of  solutions  to  models  simi¬ 
lar  to  (3.1).  For  example,  when  the  growth  rate  g  and  the  mortality  rate  m  are  functions 
of  the  size  x  only,  wellposedness  results  are  proved  in  [3]  using  weak  formulations,  and 
in  [4]  semigroup  theory  is  used  to  prove  existence  and  uniqueness  of  solutions.  For 
systems  with  time  dependent  coefficients  an  analogous  approach  using  the  theory  of 
evolution  operators  [6,  19]  offers  a  means  for  establishing  the  desired  wellposedness 
results. 


3.1  Numerical  Scheme 

The  following  notation  will  be  used  throughout  this  section:  Ax  =  xmax/nx  and  At  = 
TB/nt  denote  the  spatial  and  time  mesh  size,  respectively.  The  mesh  points  are  given 
by:  Xj  =  j Ax,  j  =  0, 1,  -  -  -  ,  nx  and  =  kAt,  k  =  0, 1,  -  -  -  ,  n*.  We  denote  by  vk , 
gk  and  mk  the  finite  difference  approximations  of  u(xj,tk),  g(xj,tk)  and  m(xj,tk), 
respectively.  We  adopt  the  following  implicit  difference  scheme  to  solve  the  model 
(3.1): 


uk+1  —  uk  ak+1uk+1  —  ak+1uk+1 

J  w  3  +  — - J  A  +  mfc+Vfc+1  =0,  1  <  j  <  nx;  0  <  A;  <  nt  —  1, 

At  Ax  J  J 

uk+1  =  0. 


(3.2) 

From  the  above  scheme  we  note  that  the  solution  remains  positive  for  all  time.  This 
is  important  from  the  biological  point  of  view. 
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3.2  Numerical  Simulations 


For  numerical  simulations  with  (3.2),  we  assume  that  xmax  =  10  gms,  Tg  =  21  days, 
i.e. ,  3  weeks.  The  mesh  size  for  Ax  and  At  are  taken  to  be  0.02  and  0.01,  respectively. 
The  initial  function  u°(x)  is  defined  as  u°(x)  =  1676505i(x).  This  implies  that  the 
initial  population  only  contains  shrimp  that  are  1  gm  in  mass,  and  the  total  number 
of  shrimp  that  are  used  to  stock  the  raceway  are 

10 

u°(x)dx  =  167650  shrimp 

which  are  1  gm  in  mass.  Hence,  the  total  biomass  is  167650  gm.  Note  that  the  stocking 
weight  of  the  shrimp  is  1  gm.  Thus,  considering  the  data  fits  of  Section  2,  we  can  assume 
that  the  growth  rate  is  a  constant  function,  and  it  is  chosen  as  g(x,t )  =  0.214.  The 
mortality  rate  function  is  chosen  as  m(x,  t )  =  0.0014.  Graphs  of  typical  solutions  are 
given  in  Figures  5  and  6.  In  Figure  5  we  plot  the  population  density  u  at  time  Tg  as  a 
function  of  the  size  x.  From  this  figure  we  see  that  at  time  t  =  Tg  the  weight  of  most 
shrimp  is  around  5.5  gm.  In  Figure  6  we  plot  the  population  density  u  as  a  function 
of  time  t  and  size  x,  where  a  cross  section  u(x,  t )  for  a  given  time  t  has  a  similar  shape 
as  u(x,Tb)  in  Figure  5  centered  at  a  different  value  of  the  size  x. 


Figure  5:  Plot  of  u(x,Tb)  as  a  function  of  the  size  x. 
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Figure  6:  Plot  of  u(x,  t )  as  a  function  of  size  x  and  time  t. 


4  Vaccine  Production  Model 

In  the  vaccine  production  stage  the  shrimp  are  infected  by  distributing  chopped  dead 
shrimp  infected  with  a  recombinant  virus  evenly  throughout  the  raceway.  This  trans¬ 
fected  biomass  is  sufficiently  large  so  that  most  of  the  shrimp  can  be  infected  in  a 
short  period  of  time,  such  as  one  day.  There  are  other  modes  of  transmission  of  virus 
in  shrimp,  such  as  cohabitation  with  infected  shrimp  that  may  be  shedding  the  virus 
into  the  surrounding  medium  (waterborne  infection).  However,  compared  to  the  prob¬ 
ability  of  shrimp  becoming  infected  via  ingestion,  these  modes  of  transmission  can  be 
assumed  (reasonably  for  this  investigation)  to  be  negligible.  Hence  we  will  only  assume 
infection  via  ingestion  of  dead  transfected  biomass.  It  is  further  assumed  that  all  the 
shrimp  have  an  equal  chance  of  becoming  infected  by  eating  the  infected  biomass.  The 
time  interval  considered  here  is  7  to  10  days.  From  [20]  and  [13]  we  know  that  during 
this  time  interval  almost  no  shrimp  progress  into  the  chronic  state.  Therefore  we  only 
consider  the  following  three  compartment  states:  susceptible  (S),  latently  infected  (L) 
and  acutely  infected  (A)  in  our  model. 

In  this  model,  we  assume  that  shrimp  will  become  instantly  infected  (i.e.,  progress 
into  latent  state)  as  soon  as  they  ingest  some  of  the  infected  biomass.  As  we  have 
noted  earlier,  however,  experimental  observations  suggest  that  there  exists  a  temporal 
delay  between  the  initial  latent  infection  and  initial  acute  infection  [14].  Moreover,  it 
is  biologically  unrealistic  to  expect  all  members  of  the  shrimp  population  to  progress 
into  the  acute  phase  at  a  fixed  number  of  days  after  initial  latent  infection.  In  addition 
the  shrimp  in  the  acute  phase  have  varying  mortality  rates  because  of  the  different 
times  that  they  progress  into  the  acute  phase  and  also  due  to  the  differences  in  genetic 
make-up  of  the  host.  As  we  have  already  noted,  we  found  that  it  is  difficult  to  account 
for  the  class  age  history  (i.e.,  the  length  of  time  that  shrimp  spend  in  a  state)  of  shrimp 
in  a  particular  (latent  or  acute)  state  using  a  system  of  delay  PDE’s  with  only  size  as 
the  structure  variable.  This  is  because  it  is  not  obvious  how  to  correctly  represent  the 
integral  involving  the  delay.  As  an  alternative,  in  order  to  model  variable  residency 
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times  we  keep  track  of  the  class  age  and  the  size  of  shrimp  by  incorporating  both  size 
structure  and  class  age  structure  into  the  latent  and  acute  states. 

We  assume  (based  on  experimental  findings)  that  there  is  a  possibility  that  shrimp 
can  stay  in  each  the  latent  and  acute  state  for  more  than  7  to  10  days.  Thus  we  can 
assume  that  the  class  age  interval  for  both  states  is  the  same  as  the  time  interval  Ty 
that  we  consider  in  our  model.  Note  that  all  shrimp  from  the  biomass  production 
raceway  are  healthy;  there  are  no  latently  infected  or  acutely  infected  shrimp  in  the 
raceway  at  time  t  =  0.  We  also  know  that  shrimp  in  the  acute  state  stop  growing, 
which  means  that  the  growth  rate  in  this  state  is  g  =  0. 

Based  on  the  above  discussions,  our  vaccine  production  model  is  given  by 

St(x,t )  +  ( gs(x)S(x,t))x  +  ms(x)S(x,t )  =  -A  S(x,t), 

Lt(x ,  t,  9)  +  ( gL(x)L(x ,  t,  6))x  +  Lg(x,  t,  9)  +  mL(x)L(x,  t,  9)  =  -7 l(9)L(x,  t,  9), 

At(x,  t,  9)  +  Ag(x,  t,  9)  +  mA(9)A(x ,  t,  9)  =  0, 

L(x,  t,  0)  =  A S(x,  t), 

A(x,t,  0)=  [  7 L(£)L(x,t,£)d£, 

Jo 

S(0,t)  =  0,  L(0,  t,  9)  =  0,  A(0,t,9)  =  0, 

S(x,  0)  =  S°(x),  L(x,0,9)  =  0,  A(x,0,6)  =  0, 

(4.1) 

where  (x,  t,  9)  G  [xmin,  x'max]  x  [0,  Ty]  x  [0,  Ty].  In  the  above  S(x,  t)  denotes  the  density 
of  individuals  (number/unit  mass)  having  mass  x  at  time  t  .  The  function  L(x,t,6) 
denotes  the  density  of  individuals  ( number /unit  mass  unit  time)  having  mass  x  at  time 
t  that  have  spent  9  days  in  the  latent  state  ,  whereas  the  function  A(x,  t,  9)  denotes 
the  density  of  individuals  (number /unit  mass  unit  time)  having  mass  x  at  time  t  that 
have  spent  9  days  in  the  acute  state  .  The  quantity  gS(x)  denotes  the  growth  rate  of 
individuals  (mass/unit  time)  in  the  susceptible  state  ,  gL(x)  denotes  the  growth  rate 
of  individuals  (mass/unit  time)  in  the  latent  state  .  The  function  ms  (x)  denotes  the 
mortality  rate  of  individuals  (1/unit  time)  in  the  susceptible  state,  and  the  function 
mL(x)  denotes  the  mortality  rate  (  1/unit  time)  of  individuals  in  the  latent  state,  and 
mA{9)  denotes  the  mortality  rate  (1/unit  time)  of  the  shrimp  that  spend  9  days  in  the 
acute  state.  The  latent  to  acute  rate  function  jL(9)  denotes  the  rate  (1/unit  time)  at 
which  the  shrimp  in  the  latent  state  that  have  spent  9  days  in  the  latent  state  become 
acutely  infected,  while  the  quantity  A  denotes  the  infection  rate  (1/unit  time)  due 
to  ingestion  of  chopped  infected  shrimp.  Finally  S°(x)  denotes  the  initial  population 
density  of  susceptible  shrimp  produced  from  the  biomass  production  model. 

4.1  Relation  to  a  Size  Structured  Epidemic  Model 

In  this  section  we  show  that  the  vaccine  production  model  (4.1)  can  be  reduced  to  a 
size  structured  three  compartment  epidemic  model. 

Let  C(x,  t )  denote  the  density  of  individuals  in  the  latent  state  having  mass  x  at 
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time  t.  Noting  that  L(x,  0,  9)  =  0,  we  have 


C(x,t)  =  f  L(x,t,6)d6.  (4.2) 

Jo 

Integrating  the  second  equation  of  (4.1)  from  0  to  t  with  respect  to  9,  we  have 

f  Lt(x,t,9)d9+  f  (gL(x)L(x,t,  0))xd0  +  f  Lg(x,t,6)d9  +  f  mL(x)L(x,t,9)d6 
Jo  Jo  Jo  Jo 

=  -  [  7 L{6)L{x,t,9)d9. 

Jo 

Simplifying  we  find 

f  Lt(x,  t,  0)dQ  +  (gL(x)C(x,  t))x  +  L(x,  t,  t )  —  L(x ,  t,  0)  +  mL(x)C(x,  t ) 

Jo 

=  -  [  'jL(9)L(x,  t,  0)d0. 

Jo 

d  fJ  L(x,  t,  9)d6  f 4 

Ct(x,t)  =  — - — - =  L(x,t,t)  +  J  Lt(x,  t,  6)d9, 

and  L(x,t,  0)  =  A S(x,t),  we  obtain 

£t(x,  t)  +  (gL(x)C(x,  t))x  +  mL{x)C(x,  t)  =  XS(x,  t)  —  f  7 L(9)L(x,t,9)d9.  (4.3) 

Jo 


If  we  define  the  average 


F  L(x,t)  = 


Jo'yL{0)L{x,t,e)d9 
J*L(x,t,6)dB  ’ 


then  on  substituting  (4.4)  in  (4.3)  we  obtain  the  equation 

£t( x,  t)  +  ( gL(x)C(x ,  t))x  +  mL(x)C(x ,  t)  =  A S(x,  t )  —  TL(x,  t)£(x,  t).  (4.5) 

Next  we  let  A(x,  t )  denote  the  density  of  individuals  in  the  acute  state  having  size 
x  at  time  t.  Noting  that  A(x,  0,  9)  =  0,  we  have 

A(x,t)  =  [  A(x,t,9)d9  (4.6) 

Jo 

Integrating  the  third  equation  of  (4.1)  from  0  to  t  with  respect  to  9,  we  have 

f  At(x,t,9)d9+  f  Ag(x,t,9)d9+  f  mA(9)A(x,  t,  9)d9  =  0. 

Jo  Jo  Jo 

Simplifying  we  have 

f  At(x,  t,  9)d9  +  A(x,  t,  t )  —  A(x,  t,  0)  +  f  mA(9)A(x,  t,  9)d9  =  0. 

Jo  Jo 


10 


Since 


and 


we  obtain 


d  A(x,t,6)d9  A 

At(x,  t)  =  — - — - =  A(x,t,t)  +  J  At(x,t,0)d0, 

A(x,  t,  0)  =  /  7  L(9)L(x,t,9)d9, 

Jo 

At(x,t)  +  [  mA(9)A(x,t,0)d9  =  f  7 L(9)L(x,t,9)d9.  (4.7) 

Jo  Jo 


>0  Jo 

As  in  the  case  for  the  latent  state,  we  define  the  average 


MA(x,t )  = 


f0  ?nA(9)A(x,  t ,  9)d9 
Jq  A(x,  t,  9)d9 


(4.8) 


Then  on  substituting  (4.8)  in  (4.7)  we  obtain  the  equation 

At(x,  t )  +  Ma(x,  t)A(x,  t )  =  rL(a:,  t)C(x,  t),  (4.9) 

Thus,  by  summing  over  the  class  age  the  vaccine  production  model  (4.1)  reduces  to 

St(x,  t )  +  ( gs(x)S(x ,  t))x  +  ms(x)S(x,  t )  =  -XS(x,  t), 

Ct(x,  t )  +  (gL(x)C(x,  t))x  +  mL(x)C(x,  t)  =  XS(x,  t )  —  TL(x,  t)C{x ,  f), 

At(.x,  t)  +  Ma(x,  t)A(x ,  t)  =  rL(a;,  t)C(x,  t ),  (4.10) 

5(0,  t)  =  0,  £(0,  t)  =  0,  A(0,  t)  =  0, 

5(x,  0)  =  5°(x),  £(x,  0)  =  0,  A(x,  0)  =  0, 

with  rL  and  M"4  as  defined  in  (4.4)  and  (4.8)  respectively.  The  model  (4.10)  is  a 
three  compartment  size  structured  model  which  involves  only  one  structure  variable, 
namely  the  mass  x. 


4.2  Method  of  Characteristics  for  the  Vaccine  Production 
Model 

In  this  section  we  use  the  method  of  characteristics  to  obtain  the  solution  of  (4.1). 
Consider  the  first  equation  of  (4.1).  The  characteristic  equations  in  this  case  are  as 
follows: 

dt{s)  -I 
ds 

t{  0)  =  t, 

which  implies  t(s)  =  t  +  s,  and 


dx(s) 

ds 


9S(x(s)), 


x(0)  =  x, 


11 


which  implies  that 


/ o  9S(x{t)) 


dx(r)  =  s.  Let  y  =  x(t),  so  that  we  have 

(  rHs)  i  \ 

1  ¥¥)iy)=eM’)- 


which  implies  that 


^Ady~ 


gs  (y) 


9s{y ) 


-s—dy  =  exp(— s). 


(4.11) 


Let  hs  (r)  =  exp 


9s (y)' 


dy  ,  then  hs (r)  is  a  decreasing  function.  Hence  its 


inverse  function  which  we  denote  by  (hb)  1  exists.  By  (4.11),  we  have  that 


hs(x(s))  =  hs(x)  exp(— s), 


which  implies  that 


x(s)  =  (hs)  1(hs(x)  exp(— s)). 


(4.12) 


We  now  consider  the  solution  of  S(x,t)  along  the  characteristic  curve.  We  then 
have 

=  _  (gr|(x(s))  +  ms(x(s))  +  A)  S(x(s),t(s)), 

S(x(0),t(0))  =  S(x,t), 
which  implies  that 


S(x,  t )  =  S(x(s),t(s))  exp  yj  ( gx  (x(t))  +  m'b(x(r))  +  A  )dr 


(4.13) 


Let  y  =  x(t)  so  that  y  =  ( hs )  1(hs(x)  exp(— r)).  We  then  find  that  r  =  In 


which  implies  that 


Noting  that  hs(y)  =  exp 


,  (hs)'(y) 

dT=~~m~v - 


c  .  .dr  ,  we  find  that 

i*-5  It)  / 


hs{y) 


(4.14) 


9S  (r) 


(F)'(y)  =  -~s 


hs{y) 
gs(y )  ‘ 


(4.15) 


By  (4.14)  and  (4.15),  we  have  dr  =  g  ^  dy.  Thus  simple  calculations  yield  that 


exp  (g% (x(r))  +  ms {x{r))dT^j  =  exp  ^ 


c/s(x(s 

— sNT  exP 
9‘ b{x) 


rx{S)  gS +  mS(y) 

»  FM  : 

I)  /  r^ms(y) 

exp  U  F(F‘ 
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Hence  from  (4.13),  we  obtain 


S(x,  t) 


9S(x(s))  ...  ,  „ 

— —  5(x(s),i(s))exp 


rx  ms(y) 
lx{a)  9S(y ) 


exp(As). 


rs  i 

Thus,  if  t  <  I  „  dy ,  we  have 

^min 


5(x,  t)  = 


gs(y) 

gs(x(-t))  co 


9s  (®) 


5u(x(— t))exp  — 


ms(y) 

'x(-t)  gsiv) 


dy  exp(— At),  (4.16) 


and  if  t  > 


gs(y ) 


dy,  we  have 


5(x,  t)  =  0. 


(4.17) 


We  then  consider  the  second  equation  of  (4.1).  As  above  we  consider  the  charac¬ 
teristic  equations 


t(s)  =  t  +  s,  9(s)  =  6  +  s 

x(s)  =  (hL)~l{hL(x)  exp(— s)),  where  hL(r)  = 


exp 


gL{y) 


dy  . 


(4.18) 


We  now  consider  the  solution  L(x,  t,  9)  along  the  characteristic  curve, 
dL(x(s),t(s),  6(s)) 


ds 


=  ~{gx(x(s))  +  m  (x(s))  +  7  {0(s))}L(x(s),t(s),9(s)) 


L(x(s),t(s),0(s))  =  L(x,t,9), 
which  implies  that 
gL(x(s )) 


L(x,  t,  9)  =  - — j-——  L(x(s),  t(s),  0(s))  exp  - 


gL(x ) 


x  mL(y ) 

(3)  gL{y) ' 


Note  that  0(s)  <  t(s)  for  every  sGt,  and  we  have  that  if  9  < 


dy,  then 


L(x,t,9 )  =  A  S(x(—9),t—9) 


gL{x{-9 )) 

gL(x) 


exp 


x  mL(y ) 

i-e)  gL{y ) 


dy  ]  exp  (  -  /  7 L{y)dyj  , 
(4.19) 


and  if  9  > 


gL{y ) 


dy,  then  we  have 


L(x,  t,  9)  =  0. 


(4.20) 


Finally,  we  consider  the  third  equation  of  (4.1).  The  characteristic  equations  in  this 
case  are 

t(s)  =  t  +  s,  9(s)  =  9  +  s,  x(s)  =  x.  (4-21) 
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We  now  consider  the  solution  A(x,  t,  9)  along  the  characteristic  curve  and  find 


dA(x(s),t(s),  0(s)) 


ds 

A(x(s),t(s),0  (s)) 


=  - mA(9(s))A(x(s),t(s),0(s )) 
A(x,t,0), 


which  implies  that 


A(x,  t,  0)  =  A(x(s),t(s),0(s))ex p  -  /  mA(y)dy  . 

V  Ms)  J 

Note  that  0(s)  <  t(s )  for  every  sGt,  and  we  have  that  if  x  >  x;mjn 

A(x,  t,  9)  =  exp  J  mA(y)dySj  J  ^L(y)L(x,t  -  9,y)dy, 


and  if  x  =  xmin,  we  have 


A(x,  t,  9)  =  0. 


(4.22) 

(4.23) 


4.3  Numerical  Scheme  for  the  Vaccine  Production  Model 

We  partition  the  x  interval  [x'min,  x'max]  with  nx  +  1  discrete  points  Xj  =  j  Ax  with 

Ax  =  m<tx - Xmm  and  j  =  0, 1, ... ,  nx.  We  divide  the  t  interval  [0,  Ty]  into  rit  +  1 

nx 

discrete  points  tk  =  kAt  with  At  =  Ty/nt  and  k  =  0,1,..., nt-  Finally  we  also 
divide  the  9  interval  [0,  Ty]  into  ng  +  1  discrete  points  9i  =  IA0  with  A 9  =  Tying  and 
Z  =  0,1,...,  ng.  We  denote  the  finite  difference  approximations  of  S(xj,tk),  L(xj,tk ,  0i) 
and  A(xj,tk,9i )  by 

S*  ~  S(xj}tk), 

Lkj’1  «  L(xj,tk,9i), 

Akj’1  «  A(xj,tk,9i), 


and  we  let 

9j  =  9S (xj),  gf  =  gL(xj ), 
nij  =  ms{xj ),  nij1  =  mL(xj) 

7^  =  7^),  mA’1  =  mA(0l). 

We  discretize  the  vaccine  production  model  (4.1)  using  an  implicit  scheme  as  follows. 
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For  j  =  1,2,...  nx,  k  =  0, 1, ...  n-t,  and  l  =  0, 1, . . . ,  n$  we  have 


(i) 

(ii) 


qk+l  _  nk  nS  qk+ 1  _  S  qk+ 1 

■  i  ,  -h  ■  i  Qj-i&j-i 


At 


+ 


Ax 


+  rrij  Sj+1  =  —XSj+1, 


-fc+i,z+i 


Y  Av-|-i y  k)l-\-\  £j  y  k~\~ lj/- 1- 1  Yj  y  k~\  l^/- )  1  t-  Aj— |— 1,Z— |— 1  t-/c-|-1,Z 

L:  ~  Li _ +  9j  Lj  -  flj-i-kj-i  fy  _  Lj 


At  Ax 

i  T  ^+4+1  _  L,Z+1  T  fe+l,i+l 

t  i 


A# 


(hi) 

(iv) 

(v) 

(Vi) 


Aj— |—  1  ,Z— |—  1  ^k,l~ (-1  Aj— |—  1  ,Z— l- 1  ^k-\-\,l 

3  3  _  _|_  3 


At 


A0  Aj -  +  mAl+1Al.+1’l+1  =  0, 


sk0+1  =  0, 

Lk0+1’l+1  =  Q)  Loi+i  =  Q)  Lk+ i,o  =  xsk+ 1; 


ne,k+ 1 


4+M+1  =  0,  A°’*+1  =  o,  Aj+1’°  =  lL’lL)+1’lAe,  where  ne,k  =  [4/A0], 

(4.24) 


z=i 


with  [z]  being  the  usual  greatest  integer  part  of  z. 

From  the  above  scheme  we  note  that  the  solution  remains  positive  for  all  time.  This 
is  important  from  the  biological  point  of  view. 


Since  C(x,t)  =  /  L(x,t,@)d0  and  A(x,t)  =  /  A(x,t,9)d0,  from  the  solution  of 
Jo  Jo 

ne,k  ne,k 

(4.24)  we  have  Ck  =  ^  Lk’lA0  and  Ak  =  ^  A1-'1  A0. 


i=i 


i=i 


Let  Bs ,  Bl  and  BA  denote  the  biomass  of  the  susceptible,  latent  and  acute  popu¬ 
lations,  respectively.  Then 


BS(t )  = 


xS(x,t)dx,  BL{t)  = 


xC(x,t)dx,  and  BA{t)  = 


xA(x,  t)dx. 


The  corresponding  ratios  of  the  biomass  in  each  state  to  the  total  biomass  are  denoted 
by 

Bs(t) 

BRatios(f)  =  BS(t)  +  BL{t)  +  BA{ty 
BRatioL(t)  =  BS(t)  +  BL{t)  +  BA{ty 

BRa,ioj(t)  =  B^t)  +  BHt)  +  BHty 

OTA 

We  also  let  N  ,  N  and  N  denote  the  total  number  in  the  susceptible,  latent  and 
acute  populations.  Then 

/#max  r^max  r^max 

S(x,t)dx,  NL(t)  =  /  £(x,t)dx,  and  NA(t)  =  /  A(x,t)dx. 

min  ”  *^min  ”  ^min 
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The  corresponding  ratios  of  the  total  number  in  each  state  to  the  total  population  are 
defined  as 

Ns(t) 

NRa*i0S<t)  =  Ns(t)  +  NL(t)  +  NA(t)' 

NL(t ) 

NRa*i0i(*)  =  Ns(t)  +  NL(t)  +  NA(t)  ’ 

NA(t) 

NRa,ioj(t)  =  wwtwwtwwy 

5  Numerical  Simulations  for  the  Coupled  Bio¬ 
mass  and  Vaccine  Model 

We  simulate  the  biomass  production  model  over  3  weeks.  The  results  of  this  model  are 
used  as  initial  conditions  to  the  vaccine  production  model.  In  the  following  examples, 
we  assume  that  xmin  =  0  grn,  xmax  =  10  grn,  Ty  =  7  days.  The  mesh  size  for  Ax, 
At  and  Ad  are  given  by  0.02,  0.01  and  0.001,  respectively.  Some  of  the  parameters  in 
(4.1)  are  defined  by  A  =  1.5,  gS(x)  =  0.214,  gL(x )  =  0.107,  ms(x)  =  0.0014,  mL(x)  = 
0.0028,  mA(0)  =  1.21  /[1.21  +  (6  —  3)2]. 

The  mortality  function  mA(0)  is  plotted  along  with  the  cumulative  mortality  in 
Figure  7.  The  cumulative  mortality  captures  the  effect  of  progression  of  TSV  infection 
on  the  mortality  rate  for  acutely  infected  shrimp.  With  this  mortality  function  we 
perform  four  different  numerical  simulations  each  with  different  choices  of  the  latent  to 
acute  rate  function.  We  present  the  resulting  plots  of  population  ratios  corresponding 
to  these  choices  of  latent  to  acute  rate  functions  to  illustrate  how  this  quantity  affects 
the  population  ratios  and  hence  the  dynamics  of  the  epidemic. 
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5.1  Example  1 

In  this  section  the  latent  to  acute  rate  function  is  taken  as 


1L(0)  =  7] i(e) 


1 

1  +  exp(-4(0  -  3)) ' 


In  Figure  8  we  plot  Ji(0)  (left)  and  the  cumulative  latent  to  acute  function  (right) 


Latent  to  Acute  Rate  Function  y)(0)  Cumulative  Function  (L->  A) 


Figure  8:  (left)  7 L{0),  (right)  1  —  exp 


7 i(T)drj.  In  this  figure  we  observe  that  after  2 

days  a  non-negligible  portion  of  latently  infected  shrimp  start  becoming  acute.  About 
2%  shrimp  stay  latently  infected  for  more  than  7  days.  In  Figure  9  we  plot  the  total 
number,  Ns ,  NL  and  NA  of  the  susceptible,  latent,  and  acute  populations,  respectively 
in  the  top  left  hand  corner.  We  plot  the  total  biomass,  BS,BL  and  BA  in  the  top 
right  hand  corner.  In  the  bottom  left  hand  corner  we  plot  the  NRatios  for  all  three 
populations,  and  in  the  bottom  right  hand  corner  we  plot  the  respective  BRatios.  All 
of  the  quantities  plotted  are  defined  in  Section  4.3.  From  this  figure  we  see  that  almost 
all  susceptible  shrimp  become  latently  infected  after  2  days.  The  maximum  number  of 
latently  infected  shrimp  is  seen  at  about  2.5  days.  We  observe  that  the  acute  phase 
starts  after  2  days.  At  day  7,  about  89%  of  shrimp  in  the  raceway  are  acutely  infected. 


which  is  given  to  be  1  —  exp  (  — 
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NRatio  Total  Number 


Figure  9:  Results  for  the  susceptible,  latent  and  acute  populations,  with 

7lw=7  m- 
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5.2  Example  2 


In  this  section  the  latent  to  acute  rate  function  is  taken  as 

.,L/q\ _ £/'" 


r(o)  =  n rw  = 


1  +  exp(-4(0  -  3))  ‘ 


In  Figure  10  we  plot  7^(0)  (left)  and  the  cumulative  latent  to  acute  function  (right) 


Latent  to  Acute  Rate  Function  7^(0) 


Cumulative  Function  (L-»  A) 


Figure  10:  (left)  7 L(0),  (right)  1  —  exp  J  7 i(r)dr 

72 ,(T)dT'Sj-  In  this  figure  we  observe 

that  around  2  days  some  portion  of  latently  infected  shrimp  start  becoming  acute. 
A  negligible  number  of  shrimp  stay  latently  infected  for  more  than  4  days.  As  done 
in  example  1,  in  Figure  11  we  plot  the  total  number,  total  biomass,  NRatios  and 
BRatios  of  the  susceptible,  latent  and  acute  populations.  We  observe  that  almost 
all  susceptible  shrimp  become  latently  infected  after  2  days.  The  maximum  number  of 
latently  infected  shrimp  occurs  at  about  2  days.  We  observe  that  the  acute  phase  starts 
around  2  days.  At  day  7,  about  99%  of  shrimp  in  the  raceway  are  acutely  infected. 


which  in  this  case  is  given  to  be  1  —  exp  I  — 


19 


Figure  11:  Results  for  the  susceptible,  latent  and  acute  populations,  with 
7LW  =  72LW- 


20 


5.3  Example  3 

In  this  section  the  latent  to  acute  rate  function  is  taken  as 


1L(0)  =  7s  (°) 


1 

1  +  exp(— 5(0  —  4)) ' 


In  Figure  12  we  plot  73  (0)  (left)  and  the  cumulative  latent  to  acute  rate  function 


Latent  to  Acute  Rate  Function  y((9) 

Cumulative  Function  (L-»  A) 


Figure  12:  (left)  7L(0),  (right)  1  —  exp 


7 |/(r)drj  .  We  observe  that  after  3  days  some  portion  of 

latently  infected  shrimp  start  becoming  acute.  About  5%  shrimp  stay  latently  infected 
for  more  than  7  days.  In  Figure  13  we  plot  the  total  number,  total  biomass,  NRatios 
and  BRatios  for  the  susceptible,  latent  and  acute  populations.  We  see  that  almost  all 
susceptible  shrimp  become  latently  infected  after  2  days.  The  maximum  number  of 
latently  infected  shrimp  occurs  at  about  3.5  days.  We  observe  that  the  acute  phase 
starts  after  3  days.  At  day  7,  about  82%  of  shrimp  in  the  raceway  are  acutely  infected. 


(right)  given  as  1  —  exp  J 
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NRatio  Total  Number 


Figure  13:  Results  for  the  susceptible,  latent,  and  acute  populations  with 
7LW  =  7aLW- 
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5.4  Example  4 

In  this  section  the  latent  to  acute  rate  function  is  taken  as 


7 L(0)  =  74  (0) 


10 

1  +  exp(— 5(0  —  4))  ’ 


In  Figure  14  we  plot  7 4  (0)  (left)  and  the  cumulative  latent  to  acute  function  (right) 


Latent  to  Acute  Rate  Function  ^(8)  Cumulative  Function  <L->  A) 


Figure  14:  (left)  7 L(0),  (right)  1  —  exp 


74  (r)rirj  .  In  this  figure  we  observe  that  around  3  days  some 

portion  of  latently  infected  shrimp  start  becoming  acute.  A  negligible  number  of  shrimp 
stay  latently  infected  for  more  than  5  days.  In  Figure  15  we  plot  the  total  number, 
total  biomass,  NRatios  and  BRatios  for  the  susceptible,  latent  and  acute  populations  as 
done  in  the  previous  examples.  Here  we  see  that  almost  all  susceptible  shrimp  become 
latently  infected  after  2  days.  The  maximum  number  of  latently  infected  shrimp  occurs 
at  about  3  days.  We  observe  that  the  acute  phase  starts  around  3  days.  At  day  7, 
about  98%  of  shrimp  in  the  raceway  are  acutely  infected. 


given  to  be  1  —  exp  — 
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NRatio  Total  Number 


Figure  15:  Results  for  the  susceptible,  latent,  and  acute  populations  with 

7LW  =  7fW- 
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5.5  Example  Summary 

In  this  section  we  compare  the  results  presented  in  examples  1,  2,  3  and  4,  to  determine 
the  effect  that  the  latent  to  acute  rate  function,  7 L(0),  has  on  the  dynamics  of  the 
epidemic.  In  Figure  16,  we  plot  the  four  different  cumulative  latent  to  acute  functions 
that  were  used  in  the  different  examples  presented  earlier.  In  Figures  17,  and  18,  we  plot 
the  total  number,  total  biomass,  NRatios  and  BRatios  for  the  simulations  performed 
with  ji,  Iz-t  and  for  the  latent  population,  and  acute  population  respectively. 

In  Figure  17  we  note  that  the  time  of  occurence  of  the  peak  of  each  curve  corre¬ 
sponding  to  total  number  or  total  biomass,  as  well  as  the  maximum  amplitude  of  each 
curve  are  determined  by  the  shape  of  the  corresponding  function  7^.  In  particular, 
the  determining  factors  are  the  time  when  a  non-negligible  number  of  latently  infected 
shrimp  start  becoming  acute,  and  the  rate  at  which  the  latent  population  becomes 
acute  as  seen  in  Figure  16. 

In  Figure  18  we  note  that  the  time  of  occurence  of  the  peak  of  each  curve  correspond¬ 
ing  to  total  number  or  total  biomass  is  determined  by  the  time  when  a  non-negligible 
number  of  latently  infected  shrimp  start  becoming  acute,  and  the  rate  at  which  the 
latent  population  becomes  acute  (see  Figure  16).  The  rate  at  which  the  population 
becomes  acute  and  the  mortality  rate  (see  Figure  7)  also  determine  the  maximum  am¬ 
plitude  of  each  curve.  For  the  acute  population  the  mortality  rate,  mA,  depends  on 
how  long  the  shrimp  stay  in  the  acute  phase  and  thus  plays  an  important  role  in  the 
dynamics  of  this  population. 

In  conclusion,  the  latent  to  acute  rate  function  7^  plays  a  very  important  role  in 
determining  the  particular  shape  of  the  curves  corresponding  to  the  total  number,  total 
biomass,  NRatio  and  BRatio,  i.e. ,  the  progression  of  the  viral  epidemic. 


Figure  16:  The  cumulative  latent  to  acute  functions. 
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NRatio  (L)  Total  Number  (L) 


Figure  17:  Results  for  the  latent  population  L. 


26 


N Ratio  (A)  Total  Number  (A) 


Figure  18:  Results  for  the  acute  population  A. 
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6  Conclusions  and  Future  Work 


In  this  paper  we  have  reported  on  our  initial  efforts  to  develop  quantitative  models 
along  with  numerical  techniques  for  simulation  of  a  coupled  biomass/ vaccine  production 
system.  While  the  overall  approach  is  valid  for  general  viruses  in  crustaceans,  we  focus 
here  on  Taura  Syndrome  Virus  in  controlled  shrimp  populations.  After  developing  the 
models  with  detailed  and  explicitly  stated  assumptions,  we  demonstrate  their  utility 
in  understanding  possible  outcomes  in  terms  of  biomass  ratios  for  several  latent  to 
acute  rate  functions  jL.  Our  numerical  simulations  in  Section  5  indicate  that  the 
distribution  of  the  infected  population  in  the  latent  and  acute  stages  is  very  sensitive 
to  the  shape  and  magnitude  of  the  function  7^.  Hence  /yL  plays  an  important  role  in 
the  construction  and  operation  of  a  production  system  such  as  the  raceways  mentioned 
in  our  introduction.  One  should  also  understand  sensitivity  of  the  model  with  respect 
to  other  functions  such  as  the  mortality  rate  mA  since  understanding  these  will  assist 
in  determining  the  optimal  harvest  time  of  infected  shrimp  in  order  to  obtain  maximal 
production  of  vaccine/ antibody.  Once  these  parameters  are  known,  one  can  use  control 
theoretic  methods  to  design  efficient  production  systems.  We  note  that,  in  general, 
parameters  such  as  7^  and  mA  will  be  species  dependent  as  well  as  virus  specific.  A 
more  formal  mathematical  sensitivity  analysis  methodology  for  the  models  we  have 
investigated  can  be  developed  using  ideas  from  [2,  5]  in  a  Prohorov  metric  framework. 

To  estimate  these  critical  parameters,  we  will  need  to  develop  an  inverse  problem 
methodology  which  will  also  help  in  ascertaining  feasibility  and  the  practical  implica¬ 
tions  of  various  operating  conditions  in  the  production  system.  We  will  also  need  data 
obtained  on  a  small  scale  from  a  population  of  shrimp  infected  with  a  specific  virus  so 
that  an  inverse  problem  can  be  performed  to  estimate  7^,  mA,  and  other  important 
parameters  of  the  production  system.  We  can  then  use  the  estimated  parameters  to 
carry  out  forward  simulations  using  the  coupled  biomass  and  vaccine  production  model 
developed  in  this  paper  to  determine  efficient  ways  to  operate  the  raceway. 

There  are  other  parameters  in  our  coupled  biomass  and  vaccine  production  model 
that  also  have  a  significant  impact  on  the  outcome  of  the  production  system.  In 
this  paper  we  have  fixed  the  growth  rates  and  the  mortality  rates  in  the  biomass 
model  as  well  as  in  all  three  compartments  in  our  vaccine  model.  This  was  done  on 
the  assumption  that  the  shrimp  populations  grow  and  are  infected  in  a  controlled 
environment.  Indeed  parameters  such  as  nutrient  levels,  temperature  and  the  density 
of  the  population,  which  can  significantly  affect  the  growth  and  mortality  rates,  can 
all  be  controlled  to  some  degree.  In  order  to  ascertain  how  these  quantities  affect  the 
dynamics  of  the  production  system  it  will  be  necessary  to  obtain  data  for  mortality 
in  shrimp  affected  by  a  particular  virus  as  well  as  growth  and  mortality  of  normal 
shrimp  in  controlled  environments.  The  infection  rate  A  may  also  be  controlled  in  the 
production  system  since  one  can  introduce  chopped  infected  biomass  into  the  raceway 
to  allow  all  the  shrimp  to  become  infected  in  a  certain  number  of  days.  In  this  paper  we 
fixed  the  infection  rate  A  to  be  a  constant  and  it  thus  is  only  determined  by  a  factor  of 
ingestion.  This  was  done  on  the  assumption  that  all  the  shrimp  have  an  equal  chance  to 
become  infected  and  all  other  modes  of  transmission  of  virus  are  negligible  compared  to 
infection  via  ingestion.  In  order  to  determine  the  optimal  amount  of  chopped  infected 
biomass  required  to  infect  the  entire  shrimp  population  in  a  certain  period  of  time,  we 
can  use  inverse  problem  techniques  along  with  data  obtained  from  experiments. 

In  the  section  on  data  fitting  we  have  shown  that  data  for  normal  growth  of  shrimp 


obtained  from  the  Waddell  Mariculture  center  could  be  fairly  approximated  with  linear 
functions  in  the  range  1-20  grn.  As  mentioned  in  Section  2,  exponential  growth  of  early 
postlarvae  has  been  observed  by  authors  in  [12].  Thus,  it  may  be  that  production  of 
vaccine  could  be  significantly  higher  from  a  given  biomass  of  postlarvae  vs.  the  same 
biomass  of  juvenile  or  adult  shrimp.  Hence  an  important  task  involves  developing  a 
second  separate  model  for  nursery  phase  of  shrimp  growth  from  zero  to  one  gram  to 
explore  use  of  very  small  animals  as  production  system  for  vaccine  or  protein.  This 
could  be  compared  to  the  existing  model.  We  would  have  to  know  the  differences  in 
dynamics  of  infection  of  small  shrimp  vs.  juveniles  and  possible  differences  in  viral 
load  per  gram  body  mass  as  shrimp  become  larger.  This  would  be  an  interesting 
application  of  the  models  for  comparison,  assuming  we  can  answer  some  of  the  basic 
questions  with  a  reporter  gene  or  endogenous  viral  protein  in  bioassay  studies.  Clearly, 
one  would  need  additional  data  in  order  to  investigate  the  importance  of  exponential 
growth  in  our  model. 

Our  model  and  associated  computational  methodology  are  sufficiently  general  to 
permit  future  investigation  of  both  underlying  crustacean  growth/death  hypotheses 
and  marine  raceway  operating  conditions  as  well  as  characteristics  of  viral  progression 
in  specific  species,  all  of  which  are  important  to  the  marine  aquaculture  industry.  For 
the  study  of  viral  epidemiology  of  marine  species,  the  calculation  of  the  basic  repro¬ 
duction  ratio,  usually  denoted  by  Ro  in  the  literature  [20],  is  important  to  determine 
whether  the  viral  disease  persists  or  dies  out.  Our  interest  in  this  paper  is  vaccine 
production,  and  calculation  of  Rq  is  not  important  in  this  context.  Computation  of  Rq 
could  in  principle,  of  course,  be  done  for  models  such  as  ours  (e.g.,  see  [16]),  but  we  do 
not  pursue  that  here. 

Investigation  of  commercial  aspects  of  the  overall  countermeasure  production  sys¬ 
tem  can  also  be  pursued  using  this  class  of  models  with  theoretic  control  methods. 
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